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BOTTOM  SEDIMENT  EFFECTS  IN  A 
NUMERICAL  PARABOLIC  EQUATION 
ENSEMBLE  AVERAGE  PROPAGATION 

Roger  M.  Oba 

Naval  Research  Laboratory, 

Bay  Saint  Louis,  MS  39529-5004  USA 

ABSTRACT 

A  method  of  numerically  computing  average  pressure  for  an  ensemble  has  been 
developed  using  a  normal  mode  analysis  of  the  parabolic  equation  [1].  It  has 
been  used  to  compute  the  transmission  loss  and  phase  for  the  average  solutions 
in  the  continuous  wave  case  where  the  vertical  variation  is  smoothly  varying.  In 
this  paper  an  alternative  choice  of  basis  functions  allows  one  to  consider  sound 
velocity  profiles  with  discontinuities  while  preserving  numerical  cc  nvergence. 
These  discontinuities  arise  when  the  sound  velocity  profile  makes  a  rapid  jump 
between  layers  of  differing  speeds.  Such  a  difference  occurs  between  the  water 
column  and  higher  velocity  fluid  bottoms  which  model  sediment  layers.  Calcu-* 
lations  using  energy  loss  due  to  sediment  penetration  allow  better  estimates  of 
attenuation. 

INTRODUCTION 

The  discussion  below  modifies  a  numerical  normal  method  for  the  propa¬ 
gation  of  the  average  of  the  full  complex  field  in  [1].  The  method  propagates 
the  complex  pressure  field  in  terms  of  a  vertical  transform  (similar  to  a  discrete 
sine  transform)  corresponding  to  eigenfunctions  of  a  piecewise  constant  sound 
velocity  profile.  First,  the  transform  is  used  to  re-express  equations  in  the  co-, 
efficient  space  of  a  reference  eigenfunction  series.  Because  this  change  bom  the 
piecewise  constant  function  to  the  refcref£e  function  involves  continuous  differen- 
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tiable  functions,  the  convergence  of  series  will  be  at  least  0(  1/n2).  Environmen¬ 
tal  fluctuation  of  sound  speed  changes  the  eigenfunctions,  and  these  changes  can 
be  expressed  as  orthogonal  matrices.  Numerical  matrix  routines  [3]  quantify  the 
changes  and  permit  their  evaluation,  as  well  as  eigenvalue  changes.  Finally,  the 
numerical  algorithm  computes  average  solutions,  taking  into  account  the  varia¬ 
tion  due  to  environmental  fluctuation. 

The  generalized  parabolic  equation  corresponds  to  the  rightward  propagat¬ 
ing  solutions  of  the  Helmholtz  equation  [2J: 

(idr  +  \Jd\  +  q  +  u  =  0.  (i; 

The  index  of  refraction  =  n  «  1,  and  the  sound  speed  term  has  jump  discon¬ 
tinuities  between  layers  of  differing  speed: 


9  =  7(2)  =  9(z)  +  <7£=o(z)  +  gt(z)  =  A:2(n2(z)-l),  (2) 

where  q(z)  is  piecewise  constant,  qt= o  and  qt  are  continuous  and  differentiable. 
For  simplicity,  the  parameter  0  <  e  <  1.  The  function  q  is  depth  dependent  and 
range  constant.  The  vertical  and  horizontal  equations  are,  respectively 


(dl  +  q+  A„)£n  =  0  (3) 

(idr  +  y/k2  -  An)  iftn  =  0.  (4) 

'  REFERENCE  VERTICAL  NORMAL  MODES 

This  section  introduces  a  vertical  sine-like  transform  of  (3)  based  on  q.  Here  we 
consider  pressure  release  top  and  bottom  surfaces.  Collins  has  already  used  sim¬ 
ilar  boundary  conditions  with  success  in  FEPE  [4].  A  single  fluid  basement  of 
higher  velocity  models  sediment,  although  several  layers  may  be  used.  (The  anal¬ 
ysis  for  a  semi-infinite  fluid  bpttom  is  nearly  identical  for  the  eigenvalues,  but 
in  the  continuous  spectrum,  it  involves  integro-differential  equations  with  Hilbert 
transforms.) 

Since  we  consider  only  two  layers  the  q  profile  is 


f  q0,  0  <  z  <  b 
\  <71,  b  <  z  <  h 


(5) 


with  q  =  qo  -  qi  >0  since  the  sediment  sound  speeds  are  higher.  There  are 
case  (a)  modes  confined  to  the  top  fluid  and  (b)  inodes  which  have  significant 


sediment  penetration. 

(a)  Top  confined  eigenvalues  have 


- 9o  <^n  <  -<7i- 


Use  fi„  =  K^  —  q 0  =  -<h  —  v\  to  define  Kn,  the  upper  layer  Form  of  the  eigenvalue, 
and  t/n,  the  lower  layer  form.  Set  z  =  h  —  z.  The  eigenfunction  and  its  derivative 


,  (  s0sinnnz,  z  <  b  '  f  s0 

^“  =  \SlstahM,  id 


SoKn  COSKnZ,  Z  <b 
-Si^n cosh l/nZ,  Z  <b 


Continuity  dictates  matching  at  layer  interface  z  =  b,  z  =  b.  Their  ratio  is 
JS-VSaS  =  -  iss^ak  .Note  that  as  v  i  0,  f  6;  so  1  <  n  <  &  +  unless 

tan  b\/q  <  -by/q  in  which  case  n  <  Determine  the  eigenvalues  numerically 
by  solving 


=  1  (■>*-' 


tanh  v„b\ 


Normalize  (where  constant  density  has  been  used)  1  =  =  So  <Pniz)dz  to 

find 

_  —\/2 
_  6sinh2»/n6  -  6  sin2  k nb  q  sin  2fcn6 

°  2sinh2i/„6  4>cnt/2 


(b)  For  sediment  penetrating  eigenvalues  q  +  n  >  0  everywhere.  Thus,  the  top 
and  bottom  forms  of  the  eigenvalues  are  defined  by  n„  =  k2  -  q0  =  -  qi . 

As  before,  the  eigenfunction  and  derivative  match  at  layer  interface  to  give  ratio 


.  Set  6  =  6  —  6.  To  find  eigenvalues  numerically,  solve 


„B  =  i  -  bl 

"  A  (1 


1  +  V  1  +  {hbq)/C 


where  £  =  nrr  4-  (-l)nSin  1  sin Kn6cosi/n6j.  This  requires  n  >  and 

n2  >  Normalization  implies  So  =  [|  +  +  a  a""2"?1"  ]  •  Finally, 

the  4>n  and  pn  depend  upon  frequency  through  the  k?  term  in  q. 

Consider  the  c  independent  case  of  q  =  q(z)  +  qo(z)  in  equation  (3)  .  It  has 
eigenfunctions  ipn(z)  and  eigenvalues  An.  Expand  these  vertical  eigenfunctions  in 
the  q  eigenfunction  scries  i/>n  =  <t>mRmn ■  This  allows  the  expression  of  the 

differential  equation  in  the  scries 

0  =  V  F  (A„  »/'))!• 

c — /m  =  1 
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Define  the  inner  product  of  velocity  and  q  eigenfunctions 
*  (gb&.tfm)  =  j°^qo{z)<f>k(*)<t>m{2)dz. 

Reformulation  of  the  original  equation  (1)  requires  definitions  of  the  several  vec¬ 
tors  and  matrices:  =  (i>k  ,  4>m )  is  the  matrix  containing  the  $  eigenfunction 

transform  of  eigenfunctions;  the  sound  velocity  profile  information  la  incorpo¬ 
rated  into  the  matrix  qo  with  component  qo.mn  =  ( )■  These  now  satisfy 
the  eigenvalue  equation 

A  =  Ar(M  -  q0)  R ,  A  =  diag(Ai,  Aa,  . . .) ,  (6) 

where  M  =  diag(pi,  pa,  . . .)  is  the  diagonal  matrix  of  q  eigenvalues. 

VERTICAL  MODES  IN  THE  ENSEMBLE  ' 

Next  take  an  ensemble  velocity  profile  expressed  in  terms  of  the  single  pa¬ 
rameter  e  in  q  =  q  +  ft  +  qt.  The  e-dependent  eigenfunctions  Vv  m  of  (3)  may  be 
written  in  V>m  series: 

,  tfre  n  =  y! .  ifimPtmni  (7) 

#m=l 

where  Ptmn  =  tym  ,  'Ptn)-  Computation  of  the  matrix  diagonalization  gives 

P[(-A  +  R\R)P,  =  -A,  (8) 


A  single  element  of  the  ensemble  (e  fixed)  is  a  solution  to 

( idr  +  y/&\  +  q  +  qo  +  uc  =  0.  (9) 

Express  the  solution  u(  of  (9)  in  an  eigenfunction  series  of  the  m  corresponding 
to  q  +  qo  +  qt  in  (3)  : 

n,  =  V  ,  0<mWem- 

"m=l 

The  coefficient  vectors  for  the  eigenfunction  series  are  related  by  RP,wt  =  u*. 
Eq.  (9)  means  vectors  corresponding  to  the  e  parameterized  eigenfunctions  sat¬ 
isfy  the  vector  differential  equation 

dPw*  =  iy/k 1  -  Atw0  (10) 

where  y/k*  —  At  is  the  diagonal  matrix  of  horizontal  eigenvalues. 
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AVERAGE  SOLUTION  TO  THE  PARABOLIC  EQUATION 

Linear  algebra  shows  that  Pf  =  eA(<)  where  A(e)  is  anti-herniitian  jfij.  Under 
suitable  conditions,  A(c)  =  cA.  In  this  ease 

P<  =  e<A  =  Ucxp(€i0)U\  (ID 

where  exp (dQ)  =  diag(e“9'  ,  ,  ...)  and  U  is  a  constant  unitary  matrix  (i 

Eq.  (10)  has  a  solution  vector 

Wt  =  exp {iyjk2  -  A«(r  -  r0])w(r0). 

This  is  orthogonally  transformed  back  to  the  q  eigenfunctions  transformed  <<>rl 
ficicnts  u(.  These  are  averaged  with  respect  to  t  to  give  the  ensemble  averaged 
with  respect  to  c 

(ut)  =  R^P,  cxp( f  \/ k2  -  A, \r  -  rol)P,T)R-Tu(r<,). 


The  parabolic  equation  (1)  has  an  ensemble  averaged  solution  expressed  in  i<  : : 
of  the  q  eigenfunction  coefficient  vectors 


(u«)  =  (RU)  (O  (RU)*u(r„), 

where 

T  =  e<'eUV'/^~^r~r,,IUf  "  * lH 
The  average  coupling  matrix  (r)  has  components 

t  OjmUjn 

with 

('y)mnj  =  J  <‘*P  ^i(0m  -  6n\  +  lyj k2  -  A,  Jr  -  r„|  j  llp{ f), 

where  dp(t)  =  '^dt  is  the  probability  density. 

ATTENUATION  ESTIMATES 

In  general,  if  the  depth  dependent  attenuation  coelicient  is  d(z),  then  the 
modal  attenuation  n~6t « lr— ro]  js  given  by 


J  u  fi 
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(Id) 


where  K,  n  =  \/k2  -  A,  [7j.  For  comparison,  f  in  [7]is  related  to  fl{z)  above  by 
inclusion  of  a  factor  of  &r  in  /3.  Computationally,  one  may  rewrite  (13)  using  (7) 
and  R  to  find 


£<n  —  -j£  ^]k kn ^j  k‘ ^ k  n 


'k]j’k 


where  Dj y  —  JQ  (3(z)(p}(z)<t>j>(z)d:.  If,  for  example,  ii  is  constant  in  the  sediment 
and  vanishes  in  the  water  column,  and  k  k'  are  in  the  non-penetrating  spec¬ 
trum,  then  Djj’  =.  /i.s-j.sy  if  we  wish  to  find 

the  ensemble  averaged  attenuation  6en, 


(&tn)  —  k'nRyk'Bj]' 

where  nfcfc-n  =  J2jy  UkjOnjU*yUny  f 
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